"""
GM1050 Advanced Principles and Algorithm of GIS
2010/2011

- Mini Project - The simplification of a map

A shapefile of departments in western Europe is given.
By the Douglas-Peuker algorithm this shapefile is simplified.
Where several topological relations should remain

@author: Bas, Josafat and Elise
"""
#-- general import
import math
import os
#from Vertex import Vertex
from dcel1 import DCEL
#from Douglas_Peucker import algorithm
from Douglas_Peucker2 import simplify_points
#-- OGR 
try:
    from osgeo import ogr
except ImportError:
    import ogr
#-- Shapely
from shapely.wkb import loads as wkbloads
from shapely.wkt import loads as wktloads

from shapely.geometry import LineString





     

def main():
    #1) Load all polygons
    #shapefile = 'final.shp';
    #shapefile = 'departementsmall.shp';
    #shapefile = 'europe_nuts.shp';
    shapefile = 'europe_nuts_mini2.shp';
    print "reading: " + shapefile
    shapelyPolygons = read_shapefile(shapefile)
    print "Number of Features: " + str(len(shapelyPolygons))
    
    #2) Extract edges and nodes - (not line segments, but a graph).
    
    #3) Structure for each polygon which edges is bounding him

    coords = createTopology(shapelyPolygons)
    
    
    text = createWkt(coords)
    outfile = raw_input('give filename (+.shp) of output file ... ')
    write_shapefile(text, outfile)
    #construct_graph(polygons)

    #total_edges = edge_counter
    #print total_edges
    #edges = graph[1]
    
    #4) Douglas - Peuker algorithm for each line
    #new_edges = {}
    #print "Give tolerance (e.g.: 0.2)"
    #tolerance = raw_input('--> ')
    #for i in range(0,total_edges):
    #    new_line = simplify_points(list(edges[i].coords), tolerance)
    #    new_edges[i] = LineString(new_line)
    
    #5) Build each polygon again, edge by edge 
    
#1) Load all polygons


def read_shapefile(infile):
    """
    Function that reads a shapefile and prints some info about it.
    
    Input:
        infile:   shapefile to read
    Output:
        - (the informations are printed in the console)

    """
    print "Reading", infile 
    ge = ogr.Open(infile)
    layer = ge.GetLayer(0) #-- a simple shapefile always has 1 layer only
    ld = layer.GetLayerDefn()
    
    #-- Find and print the geometry type of the shapefile
    print ld.GetGeomType()
    s = "Geometry type: "
    if ld.GetGeomType() == ogr.wkbPoint:
        s += "Point"
    elif ld.GetGeomType() == ogr.wkbLineString:
        s += "LineString"
    elif ld.GetGeomType() == ogr.wkbPolygon:
        s += "Polygon"
    else:
        s += "Unknown"
    print s
    
    #-- Find and print the number of geometries
    noFeat = layer.GetFeatureCount()
    print "Number of features:", noFeat
        
    #-- Find and print the attributes
    print "Attributes:"
    for i in range(ld.GetFieldCount()):
        print "\t", ld.GetFieldDefn(i).GetNameRef()
        
    #-- print the WKT of each feature
    print "The features are:"
    lsFeatures = [] #-- create an empty list of Features (of geometries)
    wktFeatures = []
    for i in xrange(noFeat):
        #print "-- --"
        f = layer.GetFeature(i)
        wkb = f.GetGeometryRef().ExportToWkb() #-- Well-Known Text
        wkt = f.GetGeometryRef().ExportToWkt() #-- Well-Known Text
        #print str(i) + "\t" + wkt 
        #wktFeatures.append({'geom': str(i) +"\t" + wkt})

        pshape = wktloads(wkt)
        
        #print str(i) + "\t" +  pshape.to_wkt()
        lsFeatures.append(pshape)
    
    return lsFeatures

#2) Extract edges and nodes - (not line segments, but a graph).

def createTopology(shapelyPolygons):
    #first we require to store information.
    #DCEL looks like the most convinient thing, but instead of storing graphs and vertexes, store nodes.
  #  print "testing"

    dcel = DCEL()
    
    dcel.convert(shapelyPolygons)
    print "----- End of conversion"
    dcel.printData();

    lineSegments = dcel.findLineSegments();
    lines = []
    tolerance = float(raw_input('Please give percentage to keep ... [0...1]'))
    #percentage = float(raw_input('Give the percentage of points to be kept [0...1] ...'))
    for lineSegmentKey in iter(lineSegments):
        lineSegment = lineSegments[lineSegmentKey]
        #print "Segment: " + lineSegmentKey 
        oldline = dcel.convertLine2Tuple(lineSegment);
        print  oldline
        newline = simplify_points(oldline, tolerance);
        print newline
        dcel.replaceLine(oldline, newline);############check for errors
        lines.append(newline)
    
    #dcel.printData();
    return lines

def createWkt(lines):
    wkt = []
    for i in range(len(lines)):
        # for each polyline
        line = ''
        for j in range(len(lines[i])):
            # for each coordinate
            line = line + str(lines[i][j][0])\
            + ' ' + str(lines[i][j][1]) + ' , '
        polyline = 'LINESTRING (' + line[0:-2] + ')'
        
        wkt.append(polyline)
    return wkt

def write_shapefile(lsFeat, outfile):
    """
    Create a new shapefile and writes it to the disk.
    
    Input:
        lsFeat:    list of Features (a dico with geom/name/pop)
        outfile:   path for the new file to create
    Output:
        - (the shapefile is written to the disk)

    """
    driver = ogr.GetDriverByName('ESRI Shapefile') #-- we create a new SHP file
    if os.path.exists(outfile):
            driver.DeleteDataSource(outfile) #-- if it exists, overwrite it
    ds = driver.CreateDataSource(outfile)
    layer = ds.CreateLayer(outfile, geom_type = ogr.wkbLineString) #-- we create a SHP with polygons
    
    #-- create the 1 attributes (called Fields in OGR jargon) "name" and "population"
    fd = ogr.FieldDefn()
    layer.CreateField(fd)
        
    for i in lsFeat:
        f = ogr.Feature(feature_def = layer.GetLayerDefn()) 
        p = ogr.CreateGeometryFromWkt(i)
        f.SetGeometry(p)
        layer.CreateFeature(f)
        f.Destroy()
    ds.Destroy()
    print "\nShapefile saved as:", outfile


def construct_graph(polygons):
    #-- several dictionaries to keep the graph organized
    # each line segment has an starting and end node
    nodes = {}
    # an linestring or linearring, describing (a part of) the boundary
    edge = {}
    # each line has one or two neighbours
    neighbour = {}
    # a check to see if the line is the boundary of an island
    island = {}
    # a counter for each edge
    global edge_counter
    edge_counter = 1
    
    #para cada poligono
    for i in range(len(polygons)):
        print "------------------------------------"
        for j in range(len(polygons)):
            if i != j:
                print str(i) + " " + str(j)
                print polygons[i]
                print polygons[j]
                intersect = polygons[i].intersection(polygons[j])
                
                test1 = polygons[i].intersection(intersect)
                test2 = polygons[j].intersection(intersect)
                print "test1: " , test1;
                print "test2: " , test2;

                #print "This polygon produces " + str(len(intersect)) + " SEGMENTS"
                #it they do not intersect, they re a geometrycollection EMPTY
                print "INTERSECT : " , intersect
                try:
                    type = intersect.geom_type
                    if type != "GeometryCollection":
                        if type != "MultiLineString":
                            print type  
                            print str(i) + " " + str(j) 
                        
                    # the common case of two neighbouring polygons result in an
                    # linearstring of linesegments
                    if type == 'GeometryCollection':
                        pass;
                    
                    if  type == 'LineString' and (((intersect.coords[0], intersect.coords[-1]) or \
                    (intersect.coords[-1], intersect.coords[0])) not in nodes):
                            #print "LineString"
                            nodes[edge_counter] = (intersect.coords[0], intersect.coords[-1])
                            edge[edge_counter] = intersect
                            neighbour[edge_counter] = (i, j)
                            island[edge_counter] = False
                            edge_counter += 1
    # a interseection can have several loose attachments to another polygon
    # this intersection results in an multi-linestring
                    elif type == 'MultiLineString':
                        pass;
                        #print "MultiLineString";
                        
#                        new_intersect = combine(intersect)
#                        for k in range(0,len(intersect)-1):
#                            if (((intersect[k].coords[0],intersect[k].coords[-1]) or \
#                            (intersect[k].coords[-1],intersect[k].coords[0])) not in nodes):
#                                nodes[edge_counter] = (intersect[k].coords[0],intersect[k].coords[-1])
#                                edge[edge_counter] = intersect[k]
#                                neighbour[edge_counter] = (i,j)
#                                island[edge_counter] = False
#                                edge_counter += 1
    # the intersection of an island with its surrounding neighbour results
    # in an linear ring
                    elif type == 'LinearRing':
                        #print "LinearRing"
                        # each island gets 3 linestrings
                        step = len(intersect.coords) / 3
                        stepping = (0, step, 2 * step, 0)
                        cor = list(intersect.coords)
                        for k in range(0, 2):
                            nodes[edge_counter] = (intersect.coords[stepping[k]], \
                            intersect.coords[stepping[k + 1]])
                            edge[edge_counter] = LineString(cor[stepping[k]:stepping[k + 1]])
                            neighbour[edge_counter] = (i, j)
                            island[edge_counter] = True
                            edge_counter += 1
                    elif type == 'MultiPolygon':
                        #print "Multipolygon"
                        print 'Not a planar partition, polygon ', i, ' & ', j, ' overlap!'
                except ValueError:
                    match = False
    #return graph



def combine(MultiLineString):
    '''
    An algorithm to combine line-segments into connected lines
    
    Input:  Shapely MultiLineString
    Output: Shapely MultiLineString
    '''
    Multistring = []
    lines = list(MultiLineString)
    for i in range(len(lines)):
        lines[i] = list(lines[i].coords)
    while len(lines) != 0:
        line = [lines[0]]
        new_line = domino(line, lines)
        lines.remove(new_line)
        Multistring.append(new_line)
    for i in len(Multistring):
        string = str(Multistring[i])
        string = string.replace("[", "(")
        string = string.replace("]", ")")
        lines.append(LineString(string))

def domino(list, rest):
    while len(rest) != 0:
        change = False
        for i in range(len(rest)):
            if list[1][1] == rest[i][0]:
                list.append(rest[i])
                rest.remove(rest[i])
                change = True
                break
            if list[0][0] == rest[i][1]:
                list.insert(0, rest[i])
                rest.remove(rest[i])
                change = True
                break
        if not change:
            break
    return list
    
#3) Structure for each polygon which edges is bounding him

#4) Douglas - Peuker algorithm for each line

#5) Build each polygon again, edge by edge 

# start of the program
if __name__ == "__main__":
    print "Running"
    main()
    #raw_input('press enter')


